Discrete rearranging disordered patterns: 
Prediction of elastic and plastic behaviour, and application to two-dimensional foams. 
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We study the elasto-plastic behaviour of materials made of individual (discrete) objects, such as 
a liquid foam made of bubbles. The evolution of positions and mutual arrangements of individual 
objects is taken into account through statistical quantities, such as the elastic strain of the structure, 
the yield strain and the yield function. The past history of the sample plays no explicit role, except 
through its effect on these statistical quantities. They suffice to relate the discrete scale with the 
collective, global scale. At this global scale, the material behaves as a continuous medium; it is 
described with tensors such as elastic strain, stress and velocity gradient. We write the differential 
equations which predict their elastic and plastic behaviour in both the general case and the case 
of simple shear. An overshoot in the shear strain or shear stress is interpreted as a rotation of the 
deformed structure, which is a purely tensorial effect that exists only if the yield strain is at least of 
order 0.3. We suggest practical applications, including: when to choose a scalar formalism rather 
than a tensorial one; how to relax trapped stresses; and how to model materials with a low, or a 
high, yield strain. 

PACS numbers: 83.80.Iz Emulsions and foams; 46.35.+Z Viscoelasticity, plasticity, viscoplasticity; 83.10.Ff 

Continuum mechanics; 47.50.-d Non-Newtonian fluid flows 

Keywords: elasticity, plasticity, overshoot, mechanics, foam, simulations 



I. INTRODUCTION 

Discrete rearranging patterns include cellular patterns, 
for instance liquid foams, biological tissues and grains in 
polycrystals; assemblies of particles such as beads, gran- 
ular materials, colloids, molecules and atoms; and inter- 
connected networks Many of these disordered ma- 
terials display elastic and plastic properties, so that the 
stress tensor can rotate and is not necessarily aligned 
with the strain rate tensor; in models this effect is in- 
cluded in objective derivatives 0. 

Use of simplified geometries, e.g. in a rhcometer, al- 
lows a first characterization of the material through mea- 
surements of shear stress. An overshoot in the shear 
stress is seen during the first loading in materials such 
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as polymers [|[ , granular materials [|| , and emulsions [|| ■ 
For liquid foams this effect has been observed in a plate- 
plate rheometer Q and in simulations 0, H| ■ It is unclear 
whether this is due to a change in the material's struc- 
ture, or a tensorial effect of shear; but nevertheless the 
overshoot is an essential ingredient in a recent model @ 
of the strain-rate discontinuity in the cylindrical Couette 
foam flow experiments of ref. (To| . Such an overshoot 
results in mechanical bistability: two different values of 
strain correspond to the same value of stress between the 
plateau and the maximum, and can thus coexist. Here, 
we investigate the elastic regime and elasto-plastic tran- 
sition in a fully tensorial model. To describe the me- 
chanical behaviour we use a formalism adapted for dis- 
crete rearranging disordered patterns which enables us 
to quantify rotational effects and to test the relevant pa- 
rameters 

We use as an example a sheared liquid foam [Hrll8| . 
Although a liquid foam consists only of gas bubbles sur- 
rounded by liquid walls, it exhibits a complex mechani- 
cal behaviour. It is elastic for small strai ns, p lastic for 
large strains and flows at large strain rates |19l - [21| . This 
behavior is useful for numerous applications such as ore 
separation, oil extraction, foods and cosmetics. The indi- 
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TABLE I. Characteristics of simulated foams. The different 
columns correspond to the symbols used in Figs. [8] and 1101 
effective liquid fraction [52], area dispersity, boundary condi- 
tions and maximal amplitude of the cycles. 

vidual objects, namely the bubbles, are easily identified, 
which makes a liquid foam (or alternatively an emulsion, 
made of droplets) a model for the study of other complex 
fluids. 

This paper is organized as follows. In Section [ill WG 
simulate the quasistatic 2D flow of a foam in a Couette 
shear geometry; we explain how we perform and rep- 
resent the measurements. In Sec. IIII1 we present our 
equations, and discuss the specific effects due to the use 
of tensors, such as the overshoot. Sec. IIVI compares 
the model and the simulation, and extracts the relevant 
information. Sec. [V] presents applications to practical 
situations, i.e. how and when to use the model. Sec. 
IVII summarizes our findings. An Appendix explains the 
notation and provides the detailed equations. 

II. SIMULATIONS 

We simulate numerically a 2D foam flowing in a lin- 
ear Couette shear geometry. Simulations of dry foams 
offer several advantages: (i) the parameters are homoge- 
neous (liquid fraction, bubble area) and controlled (no 
diffusion-driven coarsening or film rupture); (ii) the yield 
strain is of order of 0.3, which is large enough to observe 
a full tensorial elastic regime while small enough that 
plastic effects can be easily observed; (in) all physical 
quantities can be easily measured. 

A. Methods 

Several ideal, two-dimensional, dry foams (l9l | are sim- 
ulated (Fig. [Hand Table I|. We use the Surface Evolver 
[23^ in a mode in which each film is represented as a cir- 
cular arc. The value of surface tension is taken equal to 
1 throughout, without loss of generality. A realistic foam 
structure is found by minimizing the total film length 
subject to the constraint of fixed bubble areas, prescribed 
at the beginning of the simulation. The simulations are 
quasistatic, which means that the system has time to re- 
lax between successive time steps (increments in applied 
strain). Relaxation effects arc thus neglected and vis- 




FIG. 1. Example of 2D foam simulation. Pictures are suc- 
cessive snapshots of a quasi-statically sheared, fully periodic 
foam. Numbers correspond to those of Figs. [2]and[5] Bubbles 
with 6 neighbours are displayed in white, otherwise in gray. 



cosity does not need to be included. The behaviour is 
expected to be elasto-plastic. 

The simulation procedure is as follows. A Voronoi con- 
struction of randomly distributed points (2~H (not shown) 
is first used to generate a fully periodic tessellation of 
the plane. To create a confined foam, bubbles at the top 
and bottom are sequentially deleted until the required 
number of bubbles remains. In each case, the structure 
is imported into the Surface Evolver and target bub- 
ble areas prescribed, cither all the same (monodispcrse, 
5 A/ A = 0), a small random variation of up to 20% about 
monodisperse (SA/A — 0.025), or equal to the areas given 
by the Voronoi construction (SA/A = 0.66). 

The initial foam configuration for each simulation (e.g. 
label (1) in Fig. [1} is found by reducing the total film 
length to a local minimum. During this minimization, 
neighbour swappings (so-called "Tls" [l9|]) are triggered 
by deleting each film that shrinks below a certain critical 
length l c and allowing a new film to form to complete the 
process. The critical length l r defines and measures an 
effective liquid fraction, $ c ff [22| , here chosen to be very 
dry (Table [J). 
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One geometry consists of a unit cell of 400 bubbles 
with fully periodic boundary conditions to eliminate any 
artefacts due to small sample sizes. The second geome- 
try mimics more closely a real experiment, and consists 
of 296 bubbles with two parallel bars (about 15 bubble 
diameters apart) confining the foam and with periodicity 
in one direction only. 

To shear the foams, two different procedures are re- 
quired. For the periodic foams, one off-diagonal compo- 
nent of the matrix describing the periodicity of the unit 
cell is adjusted by a small amount (23|. For the confined 
foams, a small step in strain is applied by moving one of 
the confining walls a small distance and then moving all 
vertices affinely. In each case this is followed by reduction 
of the film length to a minimum correct to 16 d.p. using 
conjugate gradient (without biasing the search by intro- 
ducing any large-scale perturbations of the structure). 




1 1.5 2 
time [a.u.] 



FIG. 2. Time evolution of the reference simulation (x in Tab. 
U]). Horizontal axis: time is in arbitrary units, equivalent to 
the "cumulated strain" J \ j\dt, where t is the time and 7 
is defined up to an arbitrary prefactor; here 2.25 cycles are 
represented. Vertical axis: all curves represent U xy (left scale) 
except for the saw-tooth which is the applied 7 (right scale). 
Numbers correspond to the pictures in Fig. \T\ The first step, 
plotted with a thick solid line, starts at the □ (indicated also 
by a number 1) and its end is labelled by number 2: 7 — 
— ?> 2. The second step, plotted with a thin solid line, is 
from number 2 to 3. The third step, plotted with a middle 
solid line, starts at the (indicated also by a number 3), and 
extends to number 5: 7 = —2 — » 2. Four predictions of the 
model are plotted as dashed lines (see Fig. [5] for explanation 
of the legend); for clarity they are plotted only from 1 to 2 
and from 3 to 5: note that from 3 to 4 all predictions are 
indiscernable from the simulation. 

Each foam is subjected to at least two "saw-tooth" 
shear cycles of amplitude j m ax- Positive and negative 
steps correspond respectively to shear toward increasing 
or decreasing imposed strain 7 (Fig. [2]) 
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FIG. 3. Validation of hypotheses. Both approximations of 
section III Bl are tested during the shear of the reference sim- 
ulation (x in Fab. UJ. a) Strain eigenvalues Ui vs U%; black 
solid line: initial shear (1 — 2)which anneals the disorder; blue 
dots: shear cycles (2 — 5); dashed blue line: straight line of 
slope -1 passing through the origin, b) Deviatoric stress-strain 
relation: red, a xy vsU xy ; blue, (a xx —a yy )/2 vs (U xx — U yy )/2; 
red and blue data almost perfectly overlap (see zoom in inset); 
the slope determines 2fj,, where ^ is the shear modulus. 



B. Measurements 

At each step, the positions of the bubble centres and 
films are recorded. Tensorial quantities are measured by 
averaging over all bubbles, as follows pj. 

The texture tensor M = is computed statisti- 

cally as an average over vector links £ between centres 
of neighbouring bubbles. We assume here that the refer- 
ence texture at rest, Mo, is isotropic. We thus define it by 
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measuring the average of the determinant of M over the 
duration of the whole simulation, dot (Mo) = (det(M)). 
The elastic strain of bubbles expresses the deviation from 
the reference state, U = (logM - logM ) /2 (Eq. lAlcft . 

This tensor is symmetric by construction (U yx = U xy ); 
it can be diagonalizcd and has two eigenvalues (U\, U2) 
in two orthogonal eigendirections. The simulations (Fig. 
[3^) verify that we can reasonably assume its trace to 
be always close to zero, U\ + U2 = U xx + U yy « 0, as 
is roughly expected for an incompressible material 
Thus, due to its symmetry and vanishing trace, U has 
only two independent components: 
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(1) 



Fig. [3^, shows that our measurement of the elastic strain 
makes evident the effect of shear-induced shuffling |25|: 
the annealed foam (blue dots) differs significantly from 
the initial one (black solid line). 

The contribution to the stress of the network of bub- 
ble walls is obtained by integrating over all films [23] ; 
it yields the deviatoric {i.e. traceless) part of the elas- 
tic stress tensor a. The trace of the stress, namely the 
pressure, is unimportant here. The simulations are qua- 
sistatic and the viscous part of the stress is not relevant. 

We check (Fig. [3b) that the stress and strain are 
strongly correlated; that their correlation is linear; and 
that it is isotropic (the same for xy and xx — yy compo- 
nents) 0, [29| . Half the slope thus defines and measures 
the elastic shear modulus /z. 
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FIG. 4. Representation of the evolution of U- a) Physi- 
cal space: evolution of the point (U cos 8, U sin 8); for com- 
pleteness we also plot the opposite (and strictly equivalent) 
point (— (7 cos 8, — U sin 8). b) Component space: trajectory 
of((U xx -Uyy)/2,U xy ), that is, (U cos 29, U sin 28). 



We choose to display U only, because it is dimension- 
less, and thus more general: it makes the comparison of 
different materials easy. One possibility [H is to repre- 
sent the traceless tensor U as a circle of radius U, with 
a straight line to indicate the direction 9 of its positive 
eigenvalue: see thick lines (circle and straight line) in Fig. 
UK. We do not use it here, except in the inset of Fig. [8] 
In fact, it is easier to represent U at a given time by a 
point, enabling us to plot trajectories. Its two indepen- 
dent components can be represented in two different but 
equally useful ways, as follows. Both representations are 
equally appropriate in the problem considered here be- 
cause of the circular symmetry of the yield criterion (see 
Eq. [5}. 

First, in the case of a traceless tensor, the absolute 
value of the two eigenvalues is the same and equal to the 
amplitude U of the tensor U defined as 



U = ^/( Uxx -U yy /2) 2 + U^ (2) 

or cquivalently U = \(U 1 -U 2 )/2\ = ||U||/V2, where 

/ 2 \i/2 
||U|| = I Ejj (Uij) 1 is the euclidian norm of U. We 

call 9 the direction of the greatest eigenvalue. We call 
physical space the representation of the parameters (U,9). 
It is useful because it shows the evolution of the struc- 
ture (elongation, orientation), In particular, we plot the 
trajectory of the point ([/cos #(7), U sin 0(7)) (Fig. |4ji). 

The other representation, which has already been used 
for foams [f|, is called component space. It plots the tra- 
jectory of the point {{U xx (j) - U yy (-f))/2,U xy ('j)) (Fig. 
|4)d). It is more suitable for comparison with experimental 
data, since rheometers measure the tangential stress (xy) , 
and sometimes the normal stress difference (xx — yy). 

These two possible choices are related by 
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Complete data for one simulation are plotted in Fig. [5] 
and arc discussed in the next section. 



C. Representations 

To summarize, M, U and a characterize the current 
state of the foam. These three tensors carry here the 
same information, since ^ appears constant. In what fol- 
lows, texture, elastic strain and stress tensors are always 
aligned. 
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FIG. 5. Different representations of Fig. [2] (same symbols) according to Fig. [4] a) U xy versus 7. b) Physical space, c) U versus 
7. d) component space. The dashed circle in b) and d) has radius Uy = 0.34. 



III. MODEL 
A. Implementation 

1. Elasticity equations 

As already mentioned, we consider a quasistatic limit 
in which flow is slow enough that we may neglect viscous 
stresses. The stress is then related to the elastic strain 
(Fig [3p) . We don't consider here the effect of external 
forces, such as gravity or friction on the boundary if the 
system is confined between glass plates fl(| . 

In the present 2D case, classical plasticity [3(| sug- 
gests that the material begins to yield when the differ- 
ence between the stress eigenvalues becomes too large: 
(ci — (T2) 2 = 4oy. The yield stress ay separates a do- 
main of pure elasticity from a domain in which the mate- 
rial flows plastically. A complete set of continuous equa- 
tions (Reuss equations [30]) can then be derived; ay is 
assumed to be constant (no strain hardening). The ef- 
fect of pressure (trace of the stress) is neglected, which 



is usually a good first approximation for metals for in- 
stance [201. It must be even more appropriate for soft 
materials, like a foam, for which the shear modulus is 
several orders of magnitude smaller than the bulk modu- 
lus [l9j]. In what follows we prefer to use the component 
yield criterion [30l ] : 



(4) 



Equivalently, Eq. [4] can be written for U, since the 
deviatoric parts of a and U are proportional (Fig. [3~p). 
From Eqs. [2]and[2J we can write the yield criterion as 



U = Uy. 



(5) 



This is represented by a circle in both physical and com- 
ponent spaces (Fig. [S^d) 

For the example of foams and highly-concentrated 
emulsions, Marmottant and Graner |31( suggested that 
the transition between elastic and plastic regimes is not 
sharp, but can be described by a yield function h. This 
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function is 0, respectively 1, in the pure elastic, respec- 
tively plastic, domain. Between these two limits, both 
effects are present and the proportion seems to depend 
mostly on the clastic strain amplitude U (Fig. [B^) ■ This 
assumption was successfully tested on different flow ge- 
ometries of a 2D foam [23| . 



0.8 



0.6 



0.4 



0.2 



a)" 



0.2 0.4 0.6 0.8 

U/U Y 



0.3 
0.2 
0.1 

-0.1 
-0.2 
-0.3 



b) 



FIG. 6. (color online) Model, (a) The different yield func- 
tions h used as examples in the present paper are power laws 
h(U) = {U/UyT with n = 1 (green, dash-dots), n = 2 (red, 
dashes), n — 4 (blue, dots) and n = +oo (black, thick dashes, 
equal to everywhere except at U — Uy where it is equal 
to 1). b) Corresponding limit cycles predicted by the model, 
plotted in component space (same legend). 

By assuming that the deformation is affine [l[ and ac- 
cording to the prediction of plasticity (Eq. 22 in [29j). 
we can then write a tensorial equation of evolution of 
the texture (Eq. 9 in [l|). This dictates the evolution 
of the texture due to imposed strain (deformation, rota- 
tion) and due to relaxation (rearrangements), for details 
see Eqs. IAldlA2l 



dt 



M = M.Vv + Vv'.M 
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U 



Vv S! 



U h 



U Y J U 
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(0) 



Here d/dt is the Lagrangian derivative in time (includ- 
ing advection); M.Vv + Vv'.M is the variation of M 
due to convection by the velocity gradient Vv; as ex- 
plained in appendix IA II and IA 31 the notation U : 



is the scalar product of the elastic strain tensor 



with the symmetrized velocity gradient tensor Vv sym = 
(Vv + Vv') /2 p|; conversely, U.M is the usual prod- 



uct of tensors; here H = H 



• u 

. V 



Vv 



syrn 



is the Hcavisidc 



function, which is equal to 1 if U : Vv sym is positive and 
otherwise. 

Eq. [6j links the evolution of the foam texture with the 
clastic strain. It is quasistatic in the sense that the strain 
is relevant, not the strain-rate. It can be generalized 
to evolutions quicker than the relaxation times of the 
structure (32[. Plasticity occurs only when the elastic 
strain is oriented in the direction of shear, as expressed 
bv the Hcavisidc function H (Appendix IA 41). 



The model is continuous and analytic, without fluctu- 
ations. The information regarding disorder is recorded in 
h. Trapped stresses Q are recorded in the initial value 
Mi (or equivalcntly U,). The material's yielding crite- 
rion is encoded in Uy- The history of the material only 
plays a role in determining h, Mj (or Uj), and Uy, which 
together fully describe the material. According to the 
expression of h, Eq. [5J can be integrated analytically or 
numerically. 



2. Simple shear 

To study the structure-evolution equation, i.e. the 
competition between elasticity and plasticity, and pre- 
dict the rheological behaviour, we impose a strain rate 7 
on the material. We take x as the direction of the shear, 
which gives the following velocity field: 



Vv = 



3x $X 
3y V X 9y Vy 



= 7 




and hence 



Vv. 




(7) 



(8) 



This factor 1/2 appears when comparing the scalar and 
tensorial descriptions ( appendix IA 51) . In this geometry, 
the advection term is taken equal to zero and the result- 
ing system of equations is given in Appendix IA 41 The 
reference state Mo is considered isotropic and constant 
throughout the evolution. 

We recall that this evolution is quasistatic: 7 appears 
as a prefactor in the time evolution, Eqs. IA3I We thus 
follow the evolution with the strain 7 = J 7 dt, instead 
of the time. Tensor operations and the time evolution 
of M arc implemented by a finite difference procedure. 
Between two time steps, U is recalculated with Eq. lAlcl 



B. Predictions 

We now address the resolution of the full clasto-plastic 
set of equations IA3I Our representation underlines the 
specifically tensorial effects. 



1. Purely elastic regime 

As a first example of our representation, we consider 
here the pure elastic regime. This means that we allow 
the structure to deform elastically (stretching, contrac- 
tion), but not to relax plastically (no rearrangements). 
Our formalism allows us to describe pure elasticity by 
computing the elastic strain and its evolution when the 
material is deformed. Our formalism extends to large 
strains, even those of order one (for strains much larger 
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FIG. 7. Elastic model, represented in physical (a) and com- 
ponent (b) spaces. Different initial states are taken: Ui — 
(center, red) and Ui = 0.3 (initial points scattered around the 
circle, black). Here n +oo: all the trajectories evolve elas- 
tically. For 7 > 0, Uxy increases, that is, time evolve upwards 
on (b). If 7 < 0, these graphs are unchanged, due to their 
symmetry with respect to the horizontal axis, showing that 
the purely elastic trajectories are reversible. 
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FIG. 8. Plastic limit for 7 > 0. Model of 6 Y versus U Y (solid 
line) and corresponding representations of U as circles with 
straight lines, indicating the direction of positive eigenvalue 
(inset), for several values of Uy- Simulation points (same 
symbols as in Table [I]) are plotted for comparison. 



than one, without plasticity, the formalism of large am- 
plitude strain |33[ might be preferable). Fig. [7^b shows 
trajectories for different initial elastic strains. 

For an initially isotropic material, Ui — 0, we recover 
the classical results in the small strain limit: U ~ 7/2 
and U xy ~ 7/2. The Poynting relation [33[ thus takes 
the form of a parabola; we even extend it to an initially 
anisotropic material, U 7^ (see Fig. [T5l and |A16[) . For 
higher strains, U xy is less linear with respect to 7, due to 
the rotation of the clastic strain. 



0, h equals 1, and % equals 1: 



U 
cos 9 



= U Y , 



-4U Y 



sin# = sign(7) 
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Ve 4t ^ + 1 



(9a) 
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(9c) 



This plastic limit is represented on Fig. [8] It shows 
that the larger Uy, the less aligned U is with respect 
In all cases, U xy increases monotonically. Note that to Vv s . ym . This tensorial effect is strong because 6y de- 



this is not the case for (U xx — U yy )/2 nor for U. When 
the structure is aligned perpendicularly to the shearing 
direction (U : Vv sym < 0), it contracts (U decreases) 
under shear, until it aligns with the shear. When the 
structure is aligned parallel to the shearing direction 
(U : Vv SJ/m > 0), it stretches (U increases) under shear. 
Since U is a tensor, it can continuously decrease, change 
direction and increase again without ever vanishing (as 
opposed to a scalar, which can change sign only when it 
is equal to zero). For instance, a trajectory which starts 
with a direction opposed to that of shear has first a de- 
creasing U (contraction, with U xy negative and increas- 
ing), then an increasing U (stretching, with U xy positive 
and increasing), then a constant U (yielding, with a ro- 
tation of U towards the plastic limit). 



2. Plastic Limit 

The yield strain is the amplitude of the strain when the 
material yields, that is, a scalar number. The plastic limit 
is defined as the elastic strain tensor U obtained after an 
infinitely long shearing (7 — > +00). The amplitude of 
this tensor is that of the yield strain. Its direction is 
obtained by solving Eq. [5] when its left-hand side equals 



creases quickly with Uy. The scalar limit corresponds to 
9 ~ 45°, as discussed in Sec. IV Al 



3. Transient regime 

We now consider the shearing of an initially anisotropic 
elastic strain, Uj 7^ 0. In physical or component space, 
the state of the material is initially situated on an elastic 
trajectory and must arrive at the plastic limit point (Fig. 
[S^-d). At this limit point, an increase of elastic strain is 
immediately transformed into plastic strain. Plasticity 



may occur only if U : Vv s 



>0. 



Graphically, both in physical space and in component 
space (Fig. [§]), the plastic limit is represented by the 
point where a trajectory reaches perpendicularly the cir- 
cle at U = Uy (the elastic strain increases along the tan- 
gent of the trajectory, the plastic strain relaxes towards 
the centre of the circle: to balance each other, they must 
be parallel). 

The shape of the yield function h then determines how 
the material reaches the plastic limit. For simplicity, we 
take feasa power law function: h = (U/Uy ) n (Fig. [HJi) . 
Two examples of the resulting behaviour are plotted on 
Fig. [9] The limit n — > 00 is intuitive: the material follows 
the elastic trajectory up to U — Uy; then U is fixed and 
the plastic limit is reached by describing an arc of a circle 
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FIG. 9. Elasto-plastic model. Representation of the model 
for 7 > in physical (a,c,e) and component (b,d,f) spaces. If 
7 < 0, the vertical axis of all these graphs should be inverted, 
(a-b) For n — > +oo, all the trajectories evolve elastically (U xy 
increases, that is, time evolves in the direction of the arrows) 
up to the yield strain, here taken as Uy = 0.3, then evolve 
plastically (U constant) towards the plastic limit (blue point) 
which corresponds to a specific angle, see Eq. [9] If the tra- 
jectory reaches the yield strain on the left of the plastic limit, 
U xy passes through a maximum (overshoot), (c-d) For n — 2, 
plasticity appears more progressively, thus smoothening the 
transition between elastic and plastic regimes, and decreasing 
(or even suppressing) the overshoot, (e-f) Summary of the 
mechanical behavior: "stretching" and "contraction" accord- 
ing to the direction of U with respect to shear. Here "plas- 
ticity on" or "plasticity off" refers to the Heaviside function 
in the last term of Eq. [6] when the plasticity is progressive 
(n finite); when n increases the "plasticity on" zone narrows, 
and for n infinite it is reduced to the limit circle. 

in physical and component spaces (Fig. |^b). For other 
cases (finite n) plasticity occurs earlier (Fig. [3^f) and 
trajectories converge to the plastic limit (Fig. [9j:d). 

The behaviour changes qualitatively if the sign of 7 is 
abruptly reversed. Unlike the elastic term, the plasticity 
term is irreversible due to the Heaviside function % in 
Eq. [5] This leads to an inversion of the plastic domain 
in physical and component spaces (Fig. [Hfcf) and to a 



new plastic limit position (0y —> —0y, Eq. This 
hysteretic effect is shown on Fig. [5] by reversing 7 once 
the plastic limit is reached. For high n, the new plastic 
limit is quickly reached, since the two plastic limits are 
on the same elastic trajectory. 

If we perform alternate sign changes of 7, we observe 
that the material is stuck in a limit trajectory (Fig. [BJd) . 
This trajectory is almost insensitive to h and therefore 
close to the elastic trajectory joining the two plastic lim- 
its. This has an important consequence: once in the 
plastic regime, the elastic strain (and thus the stress) 
can not be totally relaxed if we only reverse the shearing 
direction. This is examined in more detail below (Fig. 
HU). 

4- Overshoot 

As observed for n — s- +00 (Fig. the overshoot is 
due to the transition from an clastic trajectory to the 
plastic limit. The structure itself has no overshoot: U 
increases monotonically. The overshoot appears in the 
tangential strain U xy : it is a purely tensorial effect due 
to a rotation of the structure. In fact, U xy increases in 
all elastic trajectories. Upon reaching U = Uy there is a 
sudden transition to the plastic regime. For trajectories 
to the left of the plastic limit (Fig. [3]d), we see that U xy 
decreases towards the plastic limit. The overshoot cor- 
responds to the difference between the maximum value 
of U xy (where the trajectory meets the circle) and the 
plateau value (plastic limit). 

/From Fig. [5Jd, we observe a tiny overshoot for the nor- 
mal stress difference if the trajectories reach the U = Uy 
circle to the right of the plastic limit. In that case, the 
trajectories move towards the right in the elastic regime, 
then towards the left in the plastic regime. Such tra- 
jectories correspond to the right of Fig. [9Jd, that is, a 
structure with a large trapped normal stress difference 

U XX Uyy . 

For smaller n, the elastic-plastic transition is smoother 
and the overshoot is reduced. The overshoot amplitude 
for different h and Uy is plotted on Fig. [TQ1 for the case 
of an initially isotropic structure (Ui = 0). The overshoot 
increases with Uy because the plastic limit moves away 
from the initial clastic trajectory. 

IV. COMPARISON BETWEEN SIMULATIONS 
AND THE MODEL 

A. Plastic limit 

The model fsection [III B 3|) predicts that after a few cy- 
cles a limit trajectory is reached. This trajectoryis also 
displayed in simulations by Kabla and Debregeas Q . The 
plastic limit can thus be evaluated in simulations: Uy is 
estimated by averaging U on the last plateau (using any 
of Fig. Oi-d); similarly, By is estimated by averaging 9 



9 



0.32 
0.3 

0.28 
^ 0.26 

0.24 

0.22 
0.2 



Ms 

II / 

I: / 

V. I / 
I I / 
». : / / 



0.5 



1 



1.5 

7 



2.5 



/ 'J '.S - 





f7„ 



FIG. 10. Overshoot, defined as the difference between the 
maximum of the U xy versus 7 curve, and the value of the 
plateau which follows this maximum. Calculations are per- 
formed in the case of an initially isotropic structure (Ui = 0). 
(a) Zoom over the maxima of the curves of Fig. [5^. (b) Max- 
ima for different h functions (same legend as Fig. |SJ), and 
plateau value (starred line, which is the same for all h func- 
tions), plotted vs Uy- Simulation points are plotted for com- 
parison, with the same symbols as in Table |U closed symbols 
and + correspond to the maximal value (averaged over a few 
successive points) during the first shear step; open symbols 
and x correspond to the averaged value along the plateau of 
the last shear step; Uy is measured as the plateau value of U. 
(c) Zoom of (b). 



on the last plateau of Fig. [SJd or d. On Fig. [5J results 



from simulations are compared with the model. Agree- 
ment is good, and the model captures tensorial effects, 
especially because the measured 9y deviate much from 
the 45° scalar limit. 

Taking larger l c in the simulations favours neighbour 
swappings ("Tls"), thus it corresponds to an increased 
effective liquid fraction. As expected [2(J H2 , we see a de- 
crease in Uy. However, since the effective liquid fraction 
we are simulating remains in a very dry range (< 4 10 ), 
it does not influence much Uy, which thus varies over a 
narrow range (0.26 — 0.37). 

Reaching higher Uy is possible with other materials, 
but not with disordered 2D foams. Reaching lower Uy 
is possible (and usual) in experiments on disordered wet 
foams, but not in the present simulations where the al- 
gorithm would require adaptation at high l c . 



B. Yield strain and yield function 

Simulation results fluctuate, due to the limited number 
of bubbles (discrete description) , while model curves are 
smooth, corresponding to the limit of a large number 
of bubbles (continuous material description). There is a 
qualitative agreement, which is good enough to deduce 
Uy an( l h approximately. 

For instance, Fig. [S]compares a simulation with models 
using various h functions. We observe that n ~ 2 (red 
dashes on Fig. [5]) describes well the simulation during 
the first positive shearing step, and n rts 4 (blue dotted 
lines on Fig. [5]) during the second one. Similarly, Uy is 
deduced from the plateau value of U (Figs. [51 and ITU| . 

In practice, in a first approximation, it is enough to 
consider Uy and h as constant. Their variations are small 
and thus have a small effect on the foam rhcology. How- 
ever, these variations do exist. 

For instance, in this example of Fig. [5[ n (and thus 
h) evolves throughout the simulation, revealing that the 
structure evolves too; h seems to be sensitive to the 
(topological) disorder of the foam [25|. Here Uy is con- 
stant, but there are other cases (data not shown, see pMj ) 
where, due to the decrease of the topological disorder 
during the shearing, Uy decreases. 

More generally, a real foam is constantly evolving un- 



der the effect of drainage, coarsening [21], or shuffling 
[25} . These effects should probably have to be considered 
in future models, which would try to predict Uy and h, 
based on the average and fluctuations of the structure, 
respectively. 



C. Overshoot 

We can now identify two distinct physical mechanisms 
which can cause a stress overshoot in shear experiments 
of elasto-plastic materials. 

The first one is an orientation effect, suggested in Sec. 
IIII B 41 U increases monotonically, but if Uy is large 
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enough then the rotation of U under shear implies that 
the tangential shear strain U xy passes through a maxi- 
mum. This purely tensorial effect is absent from scalar 
models. Under certain additional conditions on the ini- 
tial elastic strain Ui, which are also described correctly 
only when taking into account the tensorial aspects, the 
normal strain difference U xx — U yy too passes through a 
maximum. 

Fig. [10] shows a comparison between the model and 
the simulations. Given that in the range of simulated 
Uy the overshoot is tiny and difficult to extract from the 
fluctuations, the agreement is surprisingly good. In most 
foam experiments, where Uy is even lower, this effect 
should be too small to be measurable. 

The second one is outside of the scope of the present 
paper. It is due to an evolution of the structure itself 
during the first shear step (see section lTV B[) . This might 
be invoked to explain the larger overshoot of the data 
corresponding to the confined simulations (red and pink), 
as well as most experimental observations (such as that 
ofref. 0). 

V. PRACTICAL APPLICATIONS 

A. Comparison between scalar and tensorial 
representations 

As long as the applied shear keeps a constant direction, 
and the elastic strain remains much smaller than 1, its 
eigenvectors correspond to that of Vv sym . That is, they 
are at 45° to the direction of shear. This is called the 
scalar approximation, and it considerably simplifies the 
study of the mechanical behaviour. In that case, a single 
(scalar) number is enough to fully describe the elastic 
strain. 

This scalar number can equally well be chosen as the 
amplitude U, or the eigenvalue U±, or the tangential shear 
strain U xy , among others. To switch from one choice to 
the other requires care regarding the prcfactors [29[ : this 
is often a source of confusion in the literature, especially 
regarding the definition and value of the yield strain. The 
link between the simplified (scalar) and complete (tenso- 
rial) equations is detailed in Appendix IA 5[ using 2U xy 
as a scalar. 

If Uy <C 1, which is the case for wet foams and emul- 
sions, then U remains always much smaller than 1, and 
0y ~ 45°, so that the scalar approximation holds, see 
Fig. 19 in ref. [29| (except if the direction of the shear 
changes, in 2D or in 3D). In that limit tensorial effects 
such as normal differences or stress overshoot are negli- 
gible. 

Quantitatively, U xy is linked to sin(26>) (Eq. I3"a|) . This 
implies that a difference of 10% between the scalar and 
tensorial equations is reached when sin(2#y) = 0.9, which 
corresponds to Uy = 0.23 (Eq. [9]). Very dry foams, such 
as those simulated here, are slightly above this limit: a 
tensorial model is therefore useful. 
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2000 4000 6000 

b) time I s ] 

FIG. 11. Shearing cycles to remove trapped stresses. Here 
Ui = for the first step and Uy = 0.3, n — > +oo. Between two 
steps: the direction of shearing is turned 90° clockwise (blue 
dots); the amplitude is decreased from 7 = 2 to in 15 steps 
(black dashes); simultaneously, the direction of shear is turned 
and its amplitude is decreased (red solid line), a) Component 
space, b) (U xx — U vy )/2 versus time ( | -y | = 2.5 x 10 _3 s _1 ). 

B. Trapped strains and stresses 

A dry foam is a material with sufficiently high Uy that 
normal stresses may exist even when the material is at 
rest [H, [TH, HH. To relax such residual (or "trapped") 
stresses, we should first shear the foam enough to reach 
the plastic stage, so that plastic rearrangements anneal 
the disorder. We then must perform cycles of shear. 

If the direction of shear is kept constant, and the shear 
simply reversed, the foam asymptotically reaches a limit 
trajectory, and the stress is not relaxed. Decreasing the 
amplitude of the shear cycle does not enable to leave this 
limit trajectory (black dashes in Fig. ITTab). Kraynik 
et al. simulated dry 3D foam and applied shearing cy- 
cles (actually uniaxial contractions) of amplitude w 0.2 
in different directions, rotated by 90°; this procedure de- 
creases the trapped stress by a factor of around 2, which 
does not improve with more cycles (Fig. 7 of HH). 

Here we propose a reproducible procedure based on 
section UlI B 3[ which couples shearing cycles in different 
directions and decreasing amplitudes, as follows: 
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• The amplitude 7, of the first step is large enough 
to completely reach the plastic stage: % ^> 2Uy. 

• At each step, the shearing direction is rotated by 
90° and the shearing amplitude is decreased. 

• The decrease in amplitude between successive steps 
is smaller than 2Uy/5, ensuring there are at least 
5 steps between 2Uy and (i.e. the total number 
of steps is at least 5~fi/2Uy). 

The red solid line in Fig. ITTkb shows that the nor- 
mal stress difference decreases more at each cycle; for 
instance, 6 cycles yield a decrease by a factor of 10, ap- 
parently without saturating. In 3D the procedure is the 
same, rotating the shearing direction successively along 
the x, y and z axes [H|]. This procedure is easy to ap- 
ply to simulations, especially of fully periodic foams. In 
experiments, a special set-up should be built: in 2D, it 
can be a rubber frame in the spirit of refs. [lll.[25j. if the 
four corners can be independently displaced. 




FIG. 12. Representations of the model for Uy = 1. a) U xy 
versus 7 for the first step and different yield functions, as 
in Fig. [6ji. b) Representation of the model for 7 > in 
component spaces for n +00, as in Fig. [9b. 




FIG. 13. Representations of the model for small strains. U xy 
vs U n = (Uxx — U yy )/2. Black: exact model. Blue: first 
parabolic approximation U n = U„ + U^ y . Red: complete 
parabolic approximation (Eq. IA16[) . 



C. Materials with low and high Uy 

For practical purposes we plot the reference curves for 
two limiting types of materials: those with Uy much 
higher or much lower than 0.23. 

Fig. [T2l shows the example of Uy = 1. The plastic limit 
corresponds to a small angle 8y, resulting in a strong 
overshoot. 

Fig. [T2] shows that for small strains in the elastic 
regime, all curves can be expressed using a single pa- 
rameter; for instance, as here, which is the normal 
clastic strain U n = (U xx — U y y)/2 at zero tangential shear 
(U xy = 0). A rough parabolic approximation, and a re- 
fined one (Eq. IA16[) are plotted here for Uy = 0.3. For 
smaller Uy, this approximation is good over its whole 
range of validity (namely the elastic regime), but this 
range is smaller. 

VI. SUMMARY 

We propose a continuous model of the elasticity and 
plasticity of disordered, discrete materials such as cellu- 
lar patterns (for instance liquid foams or emulsions) and 
assemblies of particles (for instance colloids). It is based 
on statistical quantities including (i) the elastic strain U, 
a dimcnsionlcss quantity measurable on images, which 
facilitates the comparison between different experiments 
or models, and makes apparent the effect of shear on the 
material's structure; (ii) the yield strain [/, a classical 
criterion for the transition between reversible, elastic and 
irreversible, plastic regimes; (iii) and the yield function 
h(U/Uy), which describes how progressive this transi- 
tion is, by measuring the relative proportion of clastic 
and plastic deformation. They suffice to relate the dis- 
crete scale with the collective, global scale. At this global 
scale, the material behaves as a continuous medium; it is 
described with tensors such as strain, stress and velocity 
gradient. We give the differential equations which pre- 
dict the elastic and plastic behaviour. The model is fully 
tensorial and thus general, in 2D or in 3D. 

We study in detail the case of simple shear. An original 
representation, suitable for 2D incompressible materials, 
is introduced to follow the evolution of the material dur- 
ing shear. 

Since U is a tensor, it has an orientation and an am- 
plitude, which both evolve under shear. It can contin- 
uously decrease its amplitude, change direction and in- 
crease again its amplitude without ever vanishing (as op- 
posed to a scalar, which can change sign only when it 
is equal to zero). Predictions of the model regarding 
orientation and stretching are plotted. They include a 
rotation of the structure, which can induce an overshoot 
of the shear strain or shear stress (and a smaller, rarer 
overshoot in normal stress differences) even without over- 
shoot in the elastic strain amplitude. This purely tenso- 
rial effect exists if Uy is at least of order of 0.3. Indepen- 
dently, the shear can also induce a change in the mate- 
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rial's structure, sometimes resulting in a (purely scalar) 
overshoot in the modulus of the elastic strain. 

The model extends a classical plasticity criterion to dis- 
ordered media. It can be solved numerically and yields 
testable predictions. We successfully compare them with 
carefully converged quasistatic simulations of shear cycles 
in 2D foams: the clastic strain increases, saturates and 
reverses. From this comparison between model and sim- 
ulation we determine Uy and estimate h. This method 
is similar to that which we have used in experiments to 
extract Uy [3, HH, and a rough estimate of h. We still 
lack a model to predict Uy and h. Both quantities evolve 
throughout the simulation, probably due to the evolution 
of the foam's internal structure, as well as the disorder 
and fluctuations. In short, the material obeys a contin- 
uous description determined by its average properties, 
while Uy and h account for the effect at large scale of its 
fluctuations. 

All quantities involved in the model are directly mea- 
surable, as tensors, in the current state of the material; 
this includes trapped stresses which we discuss (we also 
explain how to relax them): the history of the sample 
which led to this current state plays no other direct, ex- 
plicit role. We explain how and when to use the model 
in practice, and provide a set of curves and analytical 
approximations, including a discussion and an extension 
of the Poynting relation. At low strain, typically below 
0.2, tensorial effects vanish and an approximate scalar 
simplification holds. 
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Appendix A: Detailed equations 

1. Notation for tensors 

We collect here a list of our notation, since the defini- 
tions are scattered throughout the text. For a symmetric 
tensor A, we denote by A\ and A2 its eigenvalues, by 
Ax + A 2 = A xx + A yy its trace, by A xy = A yx its off- 
diagonal term, by A n = (A xx — A yy )/2 half its normal 



difference, by A = J (A xx — A yy /2) 2 + U xy its ampli- 
tude, and by ||A|| = A\J~2 its euclidian norm defined as 
||A|| 2 = A : A. The scalar product of two tensors A and 
B is defined as A : B = J^. AijBij. 

For a traceless tensor, A\ = — A2 > 0. If 9 is the 
angle corresponding to A\, then the point with coordi- 
nates (A cos 9, A sin 9) is a representation of the actual 
direction of the tensor {physical space). Conversely, the 
point with coordinates A n = A cos 29 and A xy — A sin 29 
directly represents the components of the tensor ( compo- 
nent space). For a tensor, 9 is defined modulo 7r (and not 
2ir as for vectors), so that 29 has usually more relevance 
than 9. Similarly, for traceless tensors with eigendirec- 
tions making a relative angle <j), their scalar product is 
proportional to cos 2<j). This scalar product is maximal 
when the two eigenvectors of the positive eigenvalues co- 
incide, and minimized when they arc perpendicular. 

2. Complete system of equations 

We have obtained [H| a complete (closed) set of equa- 
tions: 

p4# = v • (-P Id + 2 / iU ) . ( Ala ) 
dt 

divu = 0, (Alb) 
U = i(bgM-IogMo), (Ale) 

— M = M.Vv + Vv'.M- 2P.M. 

dt 

(Aid) 

Eq. lAlal is the equation of dynamics, equivalent to 
Navier-Stokes, except that here the viscous stress is as- 
sumed to be negligible compared to the elastic stress. 
Eq. lAlbl assumes that the flow is incompressible; this 
assumption is often valid for foams at small deformation 
but can be relaxed if needed. Eq. lAlcl defines the elastic 
strain from the texture [l[, that is, it assumes that each 
bubble's internal degrees of freedom depend on its shape. 
Eq. lAldl is the evolution of the texture, see Eq. [5] for the 
definitions of its terms (transport and source) . Here the 
pla sticity rate P is predicted according to Eq. 22 in ref. 

p = l(w : Vv «™) n (v : Vv - m ) h (£)v- 

(A2) 

The meaning of each term is the following. The direction 
of P is set by that of U, indicating that the plasticity is 
opposed to the increase of U. The amplitude of P, that 
is the rate of plastic rearrangements, is the inverse of a 
time. It is determined by the total strain rate Vv SJ/m ; 
more precisely, by one component of Vv SJ)m , determined 
by the scalar product with U (and only if this scalar 
product is positive, as expressed by the Heaviside func- 
tion %). Finally, the amplitude of P depends on the 
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yield criterion, as expressed by H{U/Uy)' the plasticity 
appears (progressively or abruptly) when U approaches 
then exceeds the yield strain. 

Eq. IA2I is written here by assuming that M and U 
commute (see Eq. 20 of rcf. [l[), which is always the 
case if Mo is isotropic. Like Eq. IAlb| it assumes that 
the flow is incompressible, but can be extended to more 
general cases. It also assumes that the flow is slow: see 
ref. 31 1 for a discussion of "quasistatic" flow, and [HIHH 
for the extension to higher velocity. 

The next appendices examine more restrictive cases, 
that is, additional approximations: simple shear (|A 31) . 
small strain (|A 5|) . and the purely elastic regime (|A 61) . 



3. Simple shear 

In our geometry, the notation becomes: 



' sym 

M 

u 



(I 


:)■ 






M xx 


M xy 


M X y 


Myy 


U xx 


U X y\ 


U X y 


Uyy J 



U : Vv s , 



U xy j = Uj sin 26*. 



Here, due to our conventions, the angle between both 
tensors is <j> = 6 — 45°, hence the term cos 2(9 — 45°) = 
sin 28. This scalar product is maximal when the two 
eigenvectors of the positive eigenvalues coincide (which 
happens for 9 = 45°), and minimized when they are per- 
pendicular (9 = 90°). 



4. Elasto-plastic component equations 

Under simple shear the advection term is supposed 
equal to zero and Eq. lAldl becomes 



-d t M x 



2AI 



U X y 



u 



U ■ M 

- XX 

(A3a) 



d t M„ 



</</ 



U X y 

Jp' L 



U 
Uy~ 



1 

— i 

7 



-d t M xy 



M, 



yy 



H (jU xy ) h 



U ■ M 



mi 



U ■ M 



(A3b) 

J xy 

(A3c) 



The elastic regime can be studied by taking the last term 
of these equations equal to (limit of high Uy)- The 
plastic limit is calculated by taking the left hand sides of 
these equations equal to 0, h = 1, and 11=1. 



5. Scalar limit 

In the limit of small strain, U can be linearized: 
1 



U=— (M-Mo), 

ZAq 



(A4) 



where Ao is the isotropic eigenvalue of M . In that limit 
Eq. IA3cl becomes 

5) -K^) , " H0 '- ) *(£ 

(A5) 

Assuming that 9 remains close to 45° leads to 



U xx 



yy 



0. 



\u xv \ = u, 

d t 2U xy = j-jH (iUxy) h 



U x 



U Y 



(A6a) 
(A6b) 

(A6c) 



The last equation is identified as the scalar elasto-plastic 
equation [311 ] . by taking 2U xy as the scalar elastic strain. 



6. Analytical approximation at small strain 

In a purely elastic regime, the evolution equation for 
the texture is 



My 

M x 
M, 



yy 

M i -i + M i 

yy I 1 xy 

Myyl 2 + 2M xyl + Ml 



(A7) 



To express all curves analytically, we choose a single pa- 
rameter, for instance the elastic strain in a non-sheared 
state (U xy = M xy = 0): 



M % 



M x 



lvl yy I i 



x v "~yy 

which can be rewritten by eliminating 7: 



M n = Ml 



M 2 

xy 

2Mf y 



(A8) 
(A9) 

(A10) 



There are still two constants left, M l n and M yy . To elim- 
inate one of them, we use the fact that the trace of U 
is almost zero, and thus the determinant of M is almost 
constant: 



Ml Ml 



\ 2 



xx yy 

or equivalently, using Eq. IA10I 

(2M i n + Mi v )Mi n = Xl. 



(All) 



(A12) 
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Solving Eq. EH yields 




(A13) 



or equivalent ly, eliminating M yy using Eq. IA12I 
Coming back to U using Eq. IA4I 

Uxy = w^M xy , U n = -^—M n . (A15) 
ZAq ZAo 

Eq. IA14I yields a parabolic approximation: 




(A16) 



The parameter which determines each elasticity curve is 
the normal strain difference at zero shear (which is thus 
equal to the amplitude of elastic strain at zero shear) . Eq. 
IA16I is tested on Fig. [T3]for U up to 0.3. The pref actor 
of the parabola, i.e. the bracket in Eq. IA161 is exactly 1 
if XJ l n = 0: this is the Poynting relation [33[ (black curve 
on Fig. [T31 starting from the point U n = U xy = 0). In 
fact, even for XJ % n ^ 0, the bracket in Eq. IA16I remains 
close to 1: as shown in Fig. [131 the Poynting relation 
extends even to an initially anisotropic material. 
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